Issues
For overall picture look at the four plots under header PC plots
For detailed picture look at the plots under header Univariate fits and ANCOVA fits, especially compare x1 (left panel of first plot in each section) and x2 (left panel of second plot in each section)
geom_ancova_grouped <- function(m1, group_by = "pred"){
geom_smooth(method = "lm",
mapping = aes(y = predict(m1),
linetype = get(group_by))
)
}
data_folder <- "data/ghalambor"
file_name <- "allff.raw.txt"
file_path <- here(data_folder, file_name)
raw_lm <- fread(file_path) %>%
clean_names()
raw_lm[, caudal_fin := sqrt((c_fx-x6)^2 + (cfy-y6)^2)]
raw_lm[, caudal_fin2 := c_fx-x6]
# temp <- raw_lm[dataset == "greenhouse",
# .SD,
# .SDcols = c("id", "gen", "treatment", "caudal_fin", "caudal_fin2")]
raw_long <- melt(raw_lm,
id.vars = c("id", "dataset"),
measure.vars = list(paste0("x", 1:12),
paste0("y", 1:12)),
variable.name = c("landmark"),
value.name = c("x", "y")
)
qplot(x = x, y = y, data = raw_long) +
coord_fixed()
axis_levels <- c("x", "y")
lm_levels <- 1:12
lm_cols <- do.call(paste0, expand.grid(axis_levels, lm_levels))
lm_cols_x13 <- c(lm_cols, "x13")
lm_cols_13 <- c(lm_cols_x13, "y13")
lm_mat <- raw_lm[, .SD, .SDcols = lm_cols] %>%
as.matrix()
lm_array <- matrix_2_array(lm_mat)
plot_shapes(lm_array)
lm_tpr <- tpr(lm_array, 1, 6)
plot_shapes(lm_tpr)
lm_grf <- grf(lm_tpr)
plot_shapes(lm_grf)
lm_grf_dt <- array_2_matrix(lm_grf) %>%
data.table
colnames(lm_grf_dt) <- lm_cols
keep_cols <- setdiff(names(raw_lm), lm_cols)
shape <- cbind(raw_lm[, .SD, .SDcols = keep_cols],
lm_grf_dt)
setnames(shape, "treatment", "drainage_treatment")
centroid_size <- function(coords){
cols <- 1:length(coords)
even_col <- cols[lapply(cols, "%%", 2) == 0]
odd_col <- cols[lapply(cols, "%%", 2) != 0]
x <- coords[even_col]
y <- coords[odd_col]
xbar <- mean(x)
ybar <- mean(y)
cs <- sqrt(sum((x-xbar)^2 + (y-ybar)^2))
return(cs)
}
cs <- apply(lm_mat, 1, centroid_size)
shape[, centroid_size := cs]
n <- dim(lm_array)[3]
for(i in 1:n){
fig_i <- lm_array[,,i]
shape[i, median_size := median_size(fig_i)]
}
shape[, x13 := mean(x6) + caudal_fin/median_size]
data_folder <- "data/ghalambor"
file_name <- "stream to drainage map.txt"
file_path <- here(data_folder, file_name)
drainage_map <- fread(file_path) %>%
clean_names()
shape1 <- merge(shape, drainage_map,
by = "stream",
sort = FALSE)
shape <- shape1
shape[, stream5 := substr(stream_id, 1, 5) %>% tolower()]
shape[, site := paste(stream5, tolower(pred_state), sep = "-")]
# names(shape_wide)
shape_long <- melt(shape,
id.vars = c("id", "dataset", "stream_id",
"stream", "pred_state", "gen",
"drainage_treatment"),
measure.vars = list(paste0("x", 1:13),
paste0("y", 1:12)),
variable.name = c("landmark"),
value.name = c("x", "y")
)
# unique(shape_wide$landmark)
qplot(x = x, y = y, data = shape_long) +
coord_fixed()
# unique(shape$dataset)
shape[dataset == "wild", treatment := "Wild"]
shape[dataset == "garden", treatment := "No Flow"]
shape[dataset == "greenhouse" &
drainage_treatment == "pike", treatment := "Flow + pike"]
shape[dataset == "greenhouse" &
drainage_treatment == "no_pike", treatment := "Flow"]
unique(shape$stream_id)
## [1] "arimaRi" "AripoIn" "AripoDC" "ElCedroMain"
## [5] "ElCedroDavid" "ElCedroIn" "GuanapoMain" "JordanTrib"
## [9] "Marianne" "MarianneTrib" "Mauzica" "Oropuche"
## [13] "Paria" "Quare6" "Quare6_2" "QuareII"
## [17] "TurareHigh" "TurareMed" "turaretriblow" "turaremedpred"
## [21] "yaratrib" "ardc" "arin" "guan"
## [25] "orop" "quII" "turu" "quar"
shape[stream_id == "quII", stream_id := "quar"]
f2 <- shape[gen != "f0"]
# recode treatment
recode <- FALSE
if(recode == TRUE){
f2[, treatment_old := treatment]
f2[gen == "p03", treatment := "Flow + pike"]
f2[gen == "p01", treatment := "Flow"]
}
flow_levels <- c("No Flow", "Flow")
pred_levels <- c("Flow", "Flow + pike")
treatment_levels <- union(flow_levels, pred_levels)
treatment_tank_levels <- c("No Flow f2", "Flow p02", "Flow p03",
"Flow + pike p01", "Flow + pike p04")
# recoded
if(recode == TRUE){
treatment_tank_levels <- c("No Flow f2", "Flow p01", "Flow p02", "Flow + pike p03", "Flow + pike p04")
}
f2[, treatment_tank := factor(paste(treatment, gen),
levels = treatment_tank_levels)]
f2[, treatment := factor(treatment,
levels = treatment_levels)]
stream_id_levels <- c("ardc", "arin", "orop", "quar", "turu", "guan")
pred_state_levels <- c("High", "Med", "Low", "LowIntro" )
f2[, stream_id := factor(stream_id, levels = stream_id_levels)]
f2[, pred_state := factor(pred_state, levels = pred_state_levels)]
# create sample_name
f2[, sample_name := paste("fish", .I)]
# f2[, sample_name := .I]
# remove ".txt" from id
f2[, id := str_replace_all(id, ".txt", "")]
The brood file column locale is inconsistent with other files. Which is correct?
data_folder <- "data/ghalambor/rebroodvtreatment"
file_name <- "Greenhouse_Female_Data.xlsx"
file_path <- here(data_folder, file_name)
brood <- read_excel(file_path) %>%
data.table() %>%
clean_names()
# clean mark column
brood[, mark := str_replace_all(mark, " ", "")]
brood[, mark := str_replace_all(mark, ",", "_")]
brood[, mark := str_replace_all(mark, "&", "_")]
brood[, mark := str_replace_all(mark, "LLFblack", "blackLLF")]
brood[, mark := str_replace_all(mark, "and", "_")]
# what is yellowRUF_blackbelow"
# what is "purple_blackRUF" -> "purpleRUF_blackRUF"
# what is "yellowFUF" -> "yellowRUF"
# what is "greenRUF_blackLLB" -> "greenRUF_blackLLF"
brood[mark == "purple_blackRUF", mark := "purpleRUF_blackRUF"]
brood[mark == "yellowFUF", mark := "yellowRUF"]
brood[mark == "greenRUF_blackLLB", mark := "greenRUF_blackLLF"]
# Pop ID Name Locale
# 1 ardc 3
# 2 arin 4
# 3 guan 6
# 4 orop 1
# 5 quii 2
# 6 turu 5
site_list <- c("orop", "quar", "ardc", "arin", "turu", "guan")
brood[, site := site_list[locale]]
brood[, number := ifelse(id_number < 10,
paste0("0", id_number),
as.character(id_number))]
brood[, pool_id := paste0("p0", pool)]
brood[, id := paste0(site, number, pool_id)]
y_cols <- c("id", "mark")
f2 <- merge(f2, brood[, .SD, .SDcols = y_cols],
by = "id",
all.x = TRUE)
f2[, mark_all := ifelse(is.na(mark), "none", mark)]
f2[, family := tstrsplit(mark_all, "_", keep = 1)]
# f2_temp <- merge(f2, brood[, .SD, .SDcols = y_cols],
# by = "id",
# all = TRUE)
# y_cols <- c("id", "mark", "x1")
# View(f2_temp[dataset == "greenhouse", .SD, .SDcols = y_cols])
f2_marks <- f2[dataset == "greenhouse", .(N = .N),
by = .(stream_id, family, treatment_tank)]
f2_marks_wide <- dcast(f2_marks,
stream_id + family ~ treatment_tank,
value.var = "N")
PCA is not a multi-group method (at least without some modifications) so this is suboptimal but PCA should still capture most of among-group variation as long as the ratio of among:within is high.
Y <- f2[, .SD, .SDcols = lm_cols_x13] %>%
as.matrix()
row.names(Y) = f2$sample_name
S <- cov(Y)
E <- eigen(S)$vectors
f2[, pc1 := (Y %*% E[,1])[,1]]
f2[, pc2 := (Y %*% E[,2])[,1]]
f2[, pc3 := (Y %*% E[,3])[,1]]
y_cols <- c("pc1", "pc2", "pc3", lm_cols_x13)
pc_loadings <- cor(f2[, .SD, .SDcols = y_cols])[c(-1,-2,-3), 1:3]
pc_loadings %>%
kable(digits = 2) %>%
kable_styling()
| pc1 | pc2 | pc3 | |
|---|---|---|---|
| x1 | 0.48 | -0.03 | 0.11 |
| y1 | -0.03 | -0.10 | -0.07 |
| x2 | -0.79 | -0.38 | -0.28 |
| y2 | -0.12 | -0.10 | 0.34 |
| x3 | 0.45 | 0.22 | 0.29 |
| y3 | 0.23 | 0.04 | 0.68 |
| x4 | -0.81 | -0.32 | 0.26 |
| y4 | 0.32 | 0.13 | 0.51 |
| x5 | 0.39 | 0.42 | -0.33 |
| y5 | -0.48 | -0.14 | -0.20 |
| x6 | 0.44 | -0.16 | -0.56 |
| y6 | -0.56 | -0.17 | -0.48 |
| x7 | 0.12 | 0.15 | -0.74 |
| y7 | -0.17 | -0.04 | -0.50 |
| x8 | -0.36 | -0.35 | 0.31 |
| y8 | 0.04 | 0.15 | 0.05 |
| x9 | 0.39 | 0.05 | 0.44 |
| y9 | 0.34 | 0.23 | 0.07 |
| x10 | -0.03 | 0.06 | 0.31 |
| y10 | -0.04 | 0.29 | -0.18 |
| x11 | -0.06 | 0.12 | 0.09 |
| y11 | 0.39 | -0.11 | -0.07 |
| x12 | 0.10 | -0.28 | -0.10 |
| y12 | -0.33 | -0.27 | -0.09 |
| x13 | -0.86 | 0.50 | 0.00 |
gg1 <- ggplot(data = f2,
aes(x = pc1,
y = pc2,
color = treatment_tank)) +
geom_point() +
theme_pubr() +
ggtitle("PC1 v PC2: size free") +
scale_color_manual(values = pal_okabe_ito) +
theme(legend.title = element_blank()) +
guides(color = guide_legend(nrow = 2, byrow = TRUE)) +
NULL
gg2 <- ggplot(data = f2,
aes(x = pc2,
y = pc3,
color = treatment_tank)) +
geom_point() +
theme_pubr() +
ggtitle("PC2 v PC3: size free") +
scale_color_manual(values = pal_okabe_ito) +
theme(legend.title = element_blank()) +
guides(color = guide_legend(nrow = 2, byrow = TRUE)) +
NULL
plot_grid(gg1, gg2, nrow = 2)
scatter3d(pc3 ~ pc1 + pc2 | treatment_tank,
data = f2,
surface = FALSE)
# unstandardized
dist_mat <- dist(Y, method = 'euclidean')
hclust_avg <- hclust(dist_mat, method = 'average')
# plot(hclust_avg)
# extract dendrogram segment data
dendrogram_data <- dendro_data(hclust_avg)
dendrogram_segments <- dendrogram_data$segments # contains all dendrogram segment data
# get terminal dendrogram segments
dendrogram_ends <- dendrogram_segments %>%
filter(yend == 0) %>% # filter for terminal dendrogram ends
left_join(dendrogram_data$labels, by = "x") %>% # .$labels contains the row names from dist_matrix (i.e., sample_name)
rename(sample_name = label) %>%
left_join(f2, by = "sample_name")
# dataframe now contains only terminal dendrogram segments and merged metadata associated with each iris
gg_tree <- ggplot() +
geom_segment(data = dendrogram_segments,
aes(x = x,
y = y,
xend = xend,
yend = yend)) +
geom_segment(data = dendrogram_ends,
aes(x = x,
y = y.x,
xend = xend,
yend = yend,
# text = paste("sample name: ",
# sample_name, "<br>",
# "species: ", Species),
color = treatment_tank)) +
# test aes is for plotly
scale_color_manual(values = pal_okabe_ito[1:5]) +
scale_y_reverse() +
coord_flip() +
theme_bw() +
theme(legend.position = "none") +
ylab("Distance") + # flipped x and y coordinates for aesthetic reasons
ggtitle("Treatment clusters -- size free")
gg_tree
Y_raw <- f2[, .SD, .SDcols = lm_cols_x13] %>%
as.matrix()
m1 <- lm(Y ~ median_size, data = f2)
Y <- residuals(m1)
row.names(Y) = f2$sample_name
S <- cov(Y)
E <- eigen(S)$vectors
f2[, pc1_allom_free := (Y %*% E[,1])[,1]]
f2[, pc2_allom_free := (Y %*% E[,2])[,1]]
f2[, pc3_allom_free := (Y %*% E[,3])[,1]]
y_cols <- c("pc1_allom_free", "pc2_allom_free", "pc3_allom_free", lm_cols_x13)
pc_AF_loadings <- cor(f2[, .SD, .SDcols = y_cols])[c(-1,-2,-3), 1:3]
pc_both_loadings <- cbind(pc_loadings, pc_AF_loadings)
pc_both_loadings %>%
kable(digits = 2) %>%
kable_styling()
| pc1 | pc2 | pc3 | pc1_allom_free | pc2_allom_free | pc3_allom_free | |
|---|---|---|---|---|---|---|
| x1 | 0.48 | -0.03 | 0.11 | 0.18 | -0.10 | -0.32 |
| y1 | -0.03 | -0.10 | -0.07 | -0.07 | 0.07 | -0.31 |
| x2 | -0.79 | -0.38 | -0.28 | -0.68 | 0.37 | -0.31 |
| y2 | -0.12 | -0.10 | 0.34 | -0.28 | -0.07 | -0.07 |
| x3 | 0.45 | 0.22 | 0.29 | 0.40 | -0.24 | 0.24 |
| y3 | 0.23 | 0.04 | 0.68 | -0.14 | -0.29 | 0.21 |
| x4 | -0.81 | -0.32 | 0.26 | -0.72 | 0.26 | 0.42 |
| y4 | 0.32 | 0.13 | 0.51 | -0.02 | -0.34 | 0.01 |
| x5 | 0.39 | 0.42 | -0.33 | 0.53 | -0.23 | 0.05 |
| y5 | -0.48 | -0.14 | -0.20 | -0.47 | 0.08 | -0.48 |
| x6 | 0.44 | -0.16 | -0.56 | 0.50 | 0.36 | -0.14 |
| y6 | -0.56 | -0.17 | -0.48 | -0.36 | 0.25 | -0.42 |
| x7 | 0.12 | 0.15 | -0.74 | 0.28 | 0.06 | -0.39 |
| y7 | -0.17 | -0.04 | -0.50 | 0.05 | 0.19 | -0.31 |
| x8 | -0.36 | -0.35 | 0.31 | -0.42 | 0.25 | 0.24 |
| y8 | 0.04 | 0.15 | 0.05 | 0.16 | -0.05 | 0.50 |
| x9 | 0.39 | 0.05 | 0.44 | 0.15 | -0.19 | 0.15 |
| y9 | 0.34 | 0.23 | 0.07 | 0.40 | -0.12 | 0.51 |
| x10 | -0.03 | 0.06 | 0.31 | -0.06 | -0.10 | 0.41 |
| y10 | -0.04 | 0.29 | -0.18 | 0.19 | -0.13 | 0.31 |
| x11 | -0.06 | 0.12 | 0.09 | -0.01 | -0.13 | 0.07 |
| y11 | 0.39 | -0.11 | -0.07 | 0.26 | 0.13 | -0.03 |
| x12 | 0.10 | -0.28 | -0.10 | -0.08 | 0.20 | -0.39 |
| y12 | -0.33 | -0.27 | -0.09 | -0.34 | 0.23 | -0.24 |
| x13 | -0.86 | 0.50 | 0.00 | -0.61 | -0.52 | -0.10 |
cor(pc_both_loadings)[4:6, 1:3] %>%
kable(digits = 2) %>%
kable_styling()
| pc1 | pc2 | pc3 | |
|---|---|---|---|
| pc1_allom_free | 0.90 | 0.49 | -0.18 |
| pc2_allom_free | -0.30 | -0.85 | -0.50 |
| pc3_allom_free | 0.19 | 0.34 | 0.57 |
gg1 <- ggplot(data = f2,
aes(x = pc1_allom_free,
y = pc2_allom_free,
color = treatment_tank)) +
geom_point() +
theme_pubr() +
ggtitle("PC1 v PC2: allometry free") +
scale_color_manual(values = pal_okabe_ito) +
theme(legend.title = element_blank()) +
guides(color = guide_legend(nrow = 2, byrow = TRUE)) +
NULL
gg2 <- ggplot(data = f2,
aes(x = pc2,
y = pc3,
color = treatment_tank)) +
geom_point() +
theme_pubr() +
ggtitle("PC2 v PC3: allometry free") +
scale_color_manual(values = pal_okabe_ito) +
theme(legend.title = element_blank()) +
guides(color = guide_legend(nrow = 2, byrow = TRUE)) +
NULL
gg1
gg2
#plot_grid(gg1, gg2, nrow = 2)
# unstandardized
dist_mat <- dist(Y, method = 'euclidean')
hclust_avg <- hclust(dist_mat, method = 'average')
# plot(hclust_avg)
# extract dendrogram segment data
dendrogram_data <- dendro_data(hclust_avg)
dendrogram_segments <- dendrogram_data$segments # contains all dendrogram segment data
# get terminal dendrogram segments
dendrogram_ends <- dendrogram_segments %>%
filter(yend == 0) %>% # filter for terminal dendrogram ends
left_join(dendrogram_data$labels, by = "x") %>% # .$labels contains the row names from dist_matrix (i.e., sample_name)
rename(sample_name = label) %>%
left_join(f2, by = "sample_name")
# dataframe now contains only terminal dendrogram segments and merged metadata associated with each iris
gg_tree <- ggplot() +
geom_segment(data = dendrogram_segments,
aes(x = x,
y = y,
xend = xend,
yend = yend)) +
geom_segment(data = dendrogram_ends,
aes(x = x,
y = y.x,
xend = xend,
yend = yend,
# text = paste("sample name: ",
# sample_name, "<br>",
# "species: ", Species),
color = treatment_tank)) +
# test aes is for plotly
scale_color_manual(values = pal_okabe_ito[1:5]) +
scale_y_reverse() +
coord_flip() +
theme_bw() +
theme(legend.position = "none") +
ylab("Distance") + # flipped x and y coordinates for aesthetic reasons
ggtitle("Treatment clusters -- unstandardized")
gg_tree
Notes
gg1 <- ggplot(data = f2,
aes(x = median_size,
y = pc1,
color = treatment_tank)) +
geom_point() +
theme_pubr() +
ggtitle("PC1 v median size: size free") +
scale_color_manual(values = pal_okabe_ito) +
theme(legend.title = element_blank()) +
guides(color = guide_legend(nrow = 2, byrow = TRUE)) +
NULL
gg2 <- ggplot(data = f2,
aes(x = median_size,
y = pc1_allom_free,
color = treatment_tank)) +
geom_point() +
theme_pubr() +
ggtitle("PC1 v median size: allometry free") +
scale_color_manual(values = pal_okabe_ito) +
theme(legend.title = element_blank()) +
guides(color = guide_legend(nrow = 2, byrow = TRUE)) +
NULL
gg1
gg2
pd <- position_dodge(0.4)
gg <- ggplot(data = f2,
aes(x = treatment_tank,
y = pc1,
color = stream_id)) +
geom_point(position = pd) +
theme_pubr() +
scale_color_manual(values = pal_okabe_ito) +
ggtitle("Size-free PC1 by treatment") +
NULL
gg
pd <- position_dodge(0.4)
gg <- ggplot(data = f2,
aes(x = treatment_tank,
y = pc1_allom_free,
color = stream_id)) +
geom_point(position = pd) +
theme_pubr() +
scale_color_manual(values = pal_okabe_ito) +
ggtitle("Allometry-free PC1 by treatment") +
NULL
gg
pd <- position_dodge(0.4)
gg <- ggplot(data = f2,
aes(x = treatment_tank,
y = pc2,
color = stream_id)) +
geom_point(position = pd) +
theme_pubr() +
scale_color_manual(values = pal_okabe_ito) +
ggtitle("Size-free PC1 by treatment") +
NULL
gg
pd <- position_dodge(0.4)
gg <- ggplot(data = f2,
aes(x = treatment_tank,
y = pc2_allom_free,
color = stream_id)) +
geom_point(position = pd) +
theme_pubr() +
scale_color_manual(values = pal_okabe_ito) +
ggtitle("Allometry-free PC1 by treatment") +
NULL
gg
# supress message "NOTE: A nesting structure was detected in the fitted model:
# treatment_tank %in% treatment"
m1_lm <- list()
m1_lm_emm <- list()
for(j in 1:length(lm_cols_x13)){
# model fit
form_j <- paste0(lm_cols_x13[j],
" ~ stream_id * treatment_tank") %>%
formula()
m1 <- lm(form_j, data = f2)
# inference
m1_emm <- emmeans(m1,
specs = c("stream_id", "treatment_tank"))
# save to list
m1_lm[[length(m1_lm)+1]] <- m1
m1_lm_emm[[length(m1_lm_emm)+1]] <- summary(m1_emm)
}
names(m1_lm) <- lm_cols_x13
names(m1_lm_emm) <- lm_cols_x13
m1_gg <- list()
p <- length(lm_cols_x13)
pd <- position_dodge(0.4)
for(j in 1:p){
m1 <- m1_lm[[j]]
m1_emm <- m1_lm_emm[[j]]
gg <- ggplot(data = f2,
aes(x = treatment_tank,
y = get(lm_cols_x13[j]),
color = stream_id)) +
# individual points
geom_point(position = pd, alpha = 0.5) +
# group means
geom_point(data = m1_emm,
aes(x = treatment_tank,
y = emmean,
color = stream_id),
size = 3,
position = pd) +
ylab(lm_cols_x13[j]) +
theme_pubr() +
# scale_color_manual(values = pal_okabe_ito_blue) +
scale_x_discrete(labels = label_wrap(10))
NULL
m1_gg[[length(m1_gg)+1]] <- ggplotGrob(gg)
}
names(m1_gg) <- lm_cols_x13
plot_grid(m1_gg[[1]],
m1_gg[[2]], nrow = 1)
plot_grid(m1_gg[[3]],
m1_gg[[4]], nrow = 1)
plot_grid(m1_gg[[5]],
m1_gg[[6]], nrow = 1)
plot_grid(m1_gg[[7]],
m1_gg[[8]], nrow = 1)
plot_grid(m1_gg[[9]],
m1_gg[[10]], nrow = 1)
plot_grid(m1_gg[[11]],
m1_gg[[12]], nrow = 1)
plot_grid(m1_gg[[13]],
m1_gg[[14]], nrow = 1)
plot_grid(m1_gg[[15]],
m1_gg[[16]], nrow = 1)
plot_grid(m1_gg[[17]],
m1_gg[[18]], nrow = 1)
plot_grid(m1_gg[[19]],
m1_gg[[20]], nrow = 1)
plot_grid(m1_gg[[21]],
m1_gg[[22]], nrow = 1)
plot_grid(m1_gg[[23]],
m1_gg[[24]], nrow = 1)
plot_grid(m1_gg[[25]])
m1_gg <- list()
p <- length(lm_cols_x13)
pd <- position_dodge(0.4)
for(j in 1:p){
greenhouse <- f2[treatment_tank != "No Flow f2" &
stream_id != "guan",
.(y_j = mean(get(lm_cols_x13[j]))),
by = .(stream_id, treatment, treatment_tank, family)]
# model fit
# form_j <- " y_j ~ stream_id * treatment_tank + (1|family)" %>%
# formula()
# m1 <- lmer(form_j, data = greenhouse)
#
form_j <- " y_j ~ stream_id + treatment_tank + (1|family)" %>%
formula()
m1 <- lmer(form_j, data = greenhouse)
form_j <- " y_j ~ stream_id + treatment + (1|treatment_tank)" %>%
formula()
m1 <- lmer(form_j, data = greenhouse)
form_j <- " y_j ~ stream_id + treatment_tank + (treatment_tank|family)" %>%
formula()
#m1 <- aov_4(form_j, data = greenhouse)
gg <- ggplot(data = greenhouse,
aes(x = treatment_tank,
y = y_j,
color = family,
group = family)) +
# individual points
geom_point(position = pd, size = 2, show.legend = FALSE) +
ylab(lm_cols_x13[j]) +
theme_pubr() +
# scale_color_manual(values = pal_okabe_ito_blue) +
scale_x_discrete(labels = label_wrap(10))
NULL
m1_gg[[length(m1_gg)+1]] <- ggplotGrob(gg)
}
names(m1_gg) <- lm_cols_x13
n_rows <- 2
plot_grid(m1_gg[[1]],
m1_gg[[2]], nrow = n_rows)
plot_grid(m1_gg[[3]],
m1_gg[[4]], nrow = n_rows)
plot_grid(m1_gg[[5]],
m1_gg[[6]], nrow = n_rows)
plot_grid(m1_gg[[7]],
m1_gg[[8]], nrow = n_rows)
plot_grid(m1_gg[[9]],
m1_gg[[10]], nrow = n_rows)
plot_grid(m1_gg[[11]],
m1_gg[[12]], nrow = n_rows)
plot_grid(m1_gg[[13]],
m1_gg[[14]], nrow = n_rows)
plot_grid(m1_gg[[15]],
m1_gg[[16]], nrow = n_rows)
plot_grid(m1_gg[[17]],
m1_gg[[18]], nrow = n_rows)
plot_grid(m1_gg[[19]],
m1_gg[[20]], nrow = n_rows)
plot_grid(m1_gg[[21]],
m1_gg[[22]], nrow = n_rows)
plot_grid(m1_gg[[23]],
m1_gg[[24]], nrow = n_rows)
plot_grid(m1_gg[[25]])
m1_gg <- list()
p <- length(lm_cols_x13)
pd <- position_dodge(0.4)
for(j in 1:p){
greenhouse <- f2[treatment_tank != "No Flow f2",
.(y_j = mean(get(lm_cols_x13[j]))),
by = .(stream_id, treatment_tank, family)]
gg <- ggplot(data = greenhouse,
aes(x = treatment_tank,
y = y_j,
color = stream_id,
group = family)) +
# individual points
geom_point(position = pd, size = 2) +
geom_line(position = pd, alpha = 0.5) +
ylab(lm_cols_x13[j]) +
theme_pubr() +
scale_color_manual(values = pal_okabe_ito_blue) +
scale_x_discrete(labels = label_wrap(10))
NULL
m1_gg[[length(m1_gg)+1]] <- ggplotGrob(gg)
}
names(m1_gg) <- lm_cols_x13
plot_grid(m1_gg[[1]],
m1_gg[[2]], nrow = 1)
plot_grid(m1_gg[[3]],
m1_gg[[4]], nrow = 1)
plot_grid(m1_gg[[5]],
m1_gg[[6]], nrow = 1)
plot_grid(m1_gg[[7]],
m1_gg[[8]], nrow = 1)
plot_grid(m1_gg[[9]],
m1_gg[[10]], nrow = 1)
plot_grid(m1_gg[[11]],
m1_gg[[12]], nrow = 1)
plot_grid(m1_gg[[13]],
m1_gg[[14]], nrow = 1)
plot_grid(m1_gg[[15]],
m1_gg[[16]], nrow = 1)
plot_grid(m1_gg[[17]],
m1_gg[[18]], nrow = 1)
plot_grid(m1_gg[[19]],
m1_gg[[20]], nrow = 1)
plot_grid(m1_gg[[21]],
m1_gg[[22]], nrow = 1)
plot_grid(m1_gg[[23]],
m1_gg[[24]], nrow = 1)
plot_grid(m1_gg[[25]])
# supress message "NOTE: A nesting structure was detected in the fitted model:
# treatment_tank %in% treatment"
m2_lm <- list()
m2_lm_emm <- list()
for(j in 1:length(lm_cols_x13)){
# model fit
form_j <- paste0(lm_cols_x13[j],
" ~ stream_id * treatment + treatment_tank + median_size") %>%
formula()
m2 <- lm(form_j, data = f2)
# inference
# emmean is expected mean of coordinate at the average median size across all fish
m2_emm <- emmeans(m2,
specs = c("stream_id", "treatment", "treatment_tank"))
# save to list
m2_lm[[length(m2_lm)+1]] <- m2
m2_lm_emm[[length(m2_lm_emm)+1]] <- summary(m2_emm)
}
names(m2_lm) <- lm_cols_x13
names(m2_lm_emm) <- lm_cols_x13
m2_gg <- list()
p <- length(lm_cols_x13)
pd <- position_dodge(0.4)
for(j in 1:p){
m2 <- m2_lm[[j]]
m2_emm <- m2_lm_emm[[j]]
gg <- ggplot(data = f2,
aes(x = median_size,
y = get(lm_cols_x13[j]),
color = stream_id)) +
# individual points
geom_point(alpha = 0.5) +
ylab(lm_cols_x13[j]) +
theme_pubr() +
scale_color_manual(values = pal_okabe_ito_blue) +
# scale_x_discrete(labels = label_wrap(10))
NULL
m2_gg[[length(m2_gg)+1]] <- ggplotGrob(gg)
}
names(m2_gg) <- lm_cols_x13
plot_grid(m2_gg[[1]],
m2_gg[[2]], nrow = 1)
plot_grid(m2_gg[[3]],
m2_gg[[4]], nrow = 1)
plot_grid(m2_gg[[5]],
m2_gg[[6]], nrow = 1)
plot_grid(m2_gg[[7]],
m2_gg[[8]], nrow = 1)
plot_grid(m2_gg[[9]],
m2_gg[[10]], nrow = 1)
plot_grid(m2_gg[[11]],
m2_gg[[12]], nrow = 1)
plot_grid(m2_gg[[13]],
m2_gg[[14]], nrow = 1)
plot_grid(m2_gg[[15]],
m2_gg[[16]], nrow = 1)
plot_grid(m2_gg[[17]],
m2_gg[[18]], nrow = 1)
plot_grid(m2_gg[[19]],
m2_gg[[20]], nrow = 1)
plot_grid(m2_gg[[21]],
m2_gg[[22]], nrow = 1)
plot_grid(m2_gg[[23]],
m2_gg[[24]], nrow = 1)
plot_grid(m2_gg[[25]])
m2_gg <- list()
p <- length(lm_cols_x13)
pd <- position_dodge(0.4)
for(j in 1:p){
m2 <- m2_lm[[j]]
m2_emm <- m2_lm_emm[[j]]
gg <- ggplot(data = f2,
aes(x = treatment_tank,
y = get(lm_cols_x13[j]),
color = stream_id)) +
# individual points
geom_point(position = pd, alpha = 0.5) +
# group means
geom_point(data = m2_emm,
aes(x = treatment_tank,
y = emmean,
color = stream_id),
size = 3,
position = pd) +
ylab(lm_cols_x13[j]) +
theme_pubr() +
scale_color_manual(values = pal_okabe_ito_blue) +
scale_x_discrete(labels = label_wrap(10))
NULL
m2_gg[[length(m2_gg)+1]] <- ggplotGrob(gg)
}
names(m2_gg) <- lm_cols_x13
plot_grid(m2_gg[[1]],
m2_gg[[2]], nrow = 1)
plot_grid(m2_gg[[3]],
m2_gg[[4]], nrow = 1)
plot_grid(m2_gg[[5]],
m2_gg[[6]], nrow = 1)
plot_grid(m2_gg[[7]],
m2_gg[[8]], nrow = 1)
plot_grid(m2_gg[[9]],
m2_gg[[10]], nrow = 1)
plot_grid(m2_gg[[11]],
m2_gg[[12]], nrow = 1)
plot_grid(m2_gg[[13]],
m2_gg[[14]], nrow = 1)
plot_grid(m2_gg[[15]],
m2_gg[[16]], nrow = 1)
plot_grid(m2_gg[[17]],
m2_gg[[18]], nrow = 1)
plot_grid(m2_gg[[19]],
m2_gg[[20]], nrow = 1)
plot_grid(m2_gg[[21]],
m2_gg[[22]], nrow = 1)
plot_grid(m2_gg[[23]],
m2_gg[[24]], nrow = 1)
plot_grid(m2_gg[[25]])